Optimized phonon approach for the diagonalization of electron-phonon problems 
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Problems of electrons or spins interacting with lat- 
tice degrees of freedom play an important role in con- 
densed matter physics. To name only a few, consider 
for instance polaron and bi-polaron formation in various 
transition metal oxides such as tungsten oxide or high-T c 
cuprates [1], Jahn- Teller effects in colossal magnetoresis- 
tance manganites [2] , or Pcicrls and spin-Pcicrls instabil- 
ities in quasi one-dimensional materials [3]. 

As a generic model for such systems the Holstein- 
Hubbard model, 
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is frequently considered, where and b+ describe 
fermions and bosons on a site i, respectively. In many 
physically relevant situations the energy scales of all the 
subsystems - electrons (t, U), phonons (ui) and their 
interaction (gu>) - are of the same order of magnitude, 
causing analytic methods, and especially adiabatic tech- 
niques, to fail in most of these cases. Even for numerical 
methods strong interactions are a demanding task, since 
they require some cut-off in the phonon Hilbert space. 
Starting with the work of White [4] in 1993, during the 
last years a class of algorithms became very popular, 
which based on the use of a so called density matrix for 
the reduction of large Hilbert spaces to manageable di- 
mensions. Considerable focus has been placed on renor- 
malization methods for one-dimensional systems in the 
thermodynamic limit. However, exact diagonalization of 
finite clusters also benefit substantially from these ideas, 
as we will demonstrate in the present paper. 

Optimized phonon approach: First we recapitulate the 
connection between density matrices and optimized basis 
states. Starting with an arbitrary normalized quantum 
state 



M = E E ^W)\r) (2) 

r=0 u=0 

expressed in terms of the basis {|^)|r)} of the direct prod- 
uct space H = H v ® H r , we wish to reduce the dimension 
D v of the space H v by introducing a new basis, 

\v) = a ^W) > ( 3 ) 

with v = 0. . . (D v - 1) and D < D v . We call {|i>)} 
an optimized basis, if the projection of \yj) on the cor- 
responding subspace H = H„ ® H r c H is as close as 
possible to the original state. Therefore we minimize 
111-0) — \ip}\\ 2 with respect to the a^u under the condition 
(v'\v) = 5p>„, where 
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is the projected state. Since we find 
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where p = X^o" 1 llrlv'r is called the density matrix 
of the state \ip), we observe immediately that the states 
{\P)} arc optimal if they are elements of the eigenspace 
of p corresponding to its D„ largest eigenvalues wo- 

Following Zhang et al. [5], we apply these features to 
construct an optimized phonon basis for the eigenstates 
of an interacting electron/spin-phonon system. Con- 
sider a system composed of N sites, each contributing 
a phonon degree of freedom |z/j), i>i = . . . oo, and some 
other (spin or electronic) states \ri). Hence, the Hilbert 
space of the model under consideration is spanned by 
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the basis {^I^Sq 1 Wi)]^)} ■ Of course, to numerically di- 
agonalize an Hamiltonian operating on this space, we 
need to restrict ourselves to a finite-dimensional sub- 
space. To calculate, for instance, the lowest eigenstates 
of the Holstein-Hubbard model (1), we could limit the 
phonon space spanned by \vi) = (vj.)~ 1 / 2 (b\) l ' i |0) by al- 
lowing only the states < Di. Most simply we can 
choose Di — M V i yielding Z) p h = M N for the dimen- 
sion of the total phonon space. However, if we think of 
the states {(^f^ 1 as eigenstates of the Hamiltonian 
_ff p h = lo X^ilo 1 it is more suitable for most problems 
to choose an energy cut-off instead. Thus we used the 
condition Y^h=o v i < leading to D p h = ( N+ ^~ 1 ), f° r 
most of our previous numerical work (see e.g. Ref. [6]). 
For weakly interacting systems already a small num- 
ber M of phonon states is sufficient to reach very good 
convergence for ground states and low lying excitations. 
However, with increasing coupling strength most systems 
require a large number of the above 'bare' phonons, thus 
exceeding capacities of even large supercomputers. In 
some cases one can avoid these problems by choosing an 
appropriate unitary transformation of the Hamiltonian, 
but in general it is desirable to find an optimized basis 
automatically. 

The present density-matrix algorithm [5] for the con- 
struction of an optimal phonon basis considered the 
phonon subsystem as a product of one large and a 
number of small sites. Each site except the large one 
uses the same optimized basis {\fJ-i)} = {|^}} with 
v = . . . (to — 1), while the basis of the large site con- 
sists of the states plus some bare states {|^)}, 
{|/io)} = ON({|z>}} U {>}}), where ON(. . .) denotes or- 
thonormalization. After a first initialization the opti- 
mized states are improved iteratively through the follow- 
ing steps 

(1) calculating the requested eigenstate \ip) of the 
Hamiltonian H in terms of the actual basis, 

(2) replacing {l^}} with the most important (i.e. 
biggest eigenvalues wp) eigenstates of the density 
matrix p, calculated with respect to \ip) and {|/xo}}, 

(3) changing the additional states {|^)} in the set 
{|Mo}}, 

(4) orthonormalizing the set {|^o}}, and returning to 
step (1). 

A simple way to proceed in step (3) is to sweep the bare 
states {|^}} through a sufficiently large part of the infi- 
nite dimensional phonon Hilbert space. One can think of 
the algorithm as 'feeding' the optimized states with bare 
phonons, thus allowing the optimized states to become 
increasingly perfect linear combinations of bare phonon 
states. Of course the whole procedure converges only for 
eigenstates of H at the lower edge of the spectrum, since 
usually the spectrum of a Hamiltonian involving phonons 
has no upper bound. The applicability of the algorithm 



was demonstrated in Ref. [5] with the Holstein model (i.e. 
U = in Eq. (1)) as an example. 

When we implemented the above algorithm together 
with a Lanczos exact diagonalization method for our sys- 
tems of interest, we found two objections against the 
above choice of an optimized basis: (i) the basis is not 
symmetric under the symmetry operations of the Hamil- 
tonian (e.g. translations), and (ii) the phonon Hilbert 
space is still large (£> p h = Mm N ~ l , where M is the di- 
mension at the large site), since we usually need more 
than one optimized state per site. 
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FIG. 1.: Eigenvalues Wo of p calculated with the ground 
state of the Holstein model of spinless fermions for weak 
and strong coupling. For larger v the eigenvalues are 
close to the numerical precision (« 10 ), explaining 
the flattening. 

The first problem is solved by including all those states 
into the phonon basis that can be created by symmetry 
operations, and by calculating the density matrix in a 
symmetric way, i.e. by adding the density matrices gen- 
erated with respect to every site, not just site i = 0. 
Concerning the second problem we note that the eigen- 
values wo of the density matrix p decrease approximately 
exponentially, see Figure 1. If we interpret wo ~ as 
the probability of the system to occupy the corresponding 
optimized state \v), we immediately find that the prob- 
ability for the complete phonon basis state ^^Jq 1 \h) is 

proportional to £^;=o v \ This is reminiscent of the en- 
ergy cut-off discussed above, and we therefore propose 
the following choice of phonon basis states at each site, 

Vz: {| Mi )} = ON({| /U )}) (6) 
. . f opt. state \u), < n < to . . 

^' ~ { bare state \v), m < fi < M ' ^ ' 

and for the complete phonon basis |(8)s iMi <Af Im»)}> 

yielding _D p h = ( N+ ^f 1 ). Implementation of this op- 
timization procedure together with our existing Lanczos 
diagonalization code [6] allows the study of interacting 
electron/spin-phonon systems in a much larger parame- 
ter space without reaching the limits of available super- 
computers. 
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To demonstrate the power of the method, in the fol- 
lowing we adress two frequently discussed problems: bi- 
polaron formation in the Holstein-Hubbard model, and 
Luttinger liquid characteristics of a ID polaronic metal. 
(a) Bi-polarons in the Holstein-Hubbard model have been 
the subject of numerous studies over the last decades, 
stimulated for instance by the discovery of high-T c 
cuprates, and the belief that the interplay between strong 
electron-phonon and electron-electron interactions plays 
a significant role in these highly correlated materials [7] . 
Nevertheless the influence of the Hubbard interaction U 
on bi-polaron formation is still not completely under- 
stood. Beside bi-polaron formation itself, an interesting 
open question is the transition between two bi-polaronic 
regimes, namely the inter-site and the on-site bi-polaron. 
Since the Hubbard interaction U and the electron-phonon 
interaction compete, we usually need to consider inter- 
mediate to strong electron-phonon coupling g, or A := 
g 2 u/(2t), making the problem a good testing ground for 
our optimized phonon algorithm. 
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FIG. 3.: Bi-polaron kinetic energy as a function of cou- 
pling strength A for frequencies Lu/t = 0.4 and 3.0, com- 
paring optimized (OPM) and bare (ED) phonons. Inset: 
bi-polaron dispersion at u/t = 0.4, A = U/(4t), where 
AE/t = 0.0103 is the bandwidth. 
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FIG. 2.: Critical coupling for inter-site bi-polaron forma- 
tion at fixed U/t — 6.0 and system size N — 12; inset: 
(niUj) correlation in the unbound, inter-site and on-site 
cases for u/t = 0.4, N = 6. 

In a recent work Bonca et al. [8] studied mobile bi- 
polarons in the Hubbard model with the aid of a varia- 
tional technique. Their focus is mainly on the U depen- 
dence of the transition from unbound polarons to inter- 
site bi-polarons and from inter-site to on-site bi-polarons 
at intermediate frequencies. These transitions also show 
a significant u dependence. Here the adiabatic frequency 
range is of special interest, since there are no appropriate 
analytic methods for small but finite frequencies. 

Using the optimized phonon approach on lattices sizes 
up to iV = 12, we calculated the phase diagram for 
the transition from unbound polarons to inter-site bi- 
polarons at fixed U. The critical coupling A c was de- 
termined by the condition A = 0, where A is the en- 
ergy difference between the two-particle ground state and 
twice the one-particle ground state, i.e. A = Eb — 2E p . 
As indicated in Figure 2, the critical interaction A c in- 



creases with frequency, reaching U/{At) as the limiting 
value. The density-density correlation (see inset) signals 
a second transition from inter-site to on-site bi-polarons, 
which also causes a distinct feature in the kinetic en- 
ergy [9], shown in Figure 3 for different frequencies and 
system sizes. For u — OAt there is a sharp drop in £kin 
at about A = 1.6. Near this critical coupling we ob- 
serve another striking effect if we study the bi-polaron 
band dispersion: Namely, it is almost a perfect cosine at 
the critical coupling, but deviates from this simple tight- 
binding dispersion for other couplings. That means in 
the vicinty of A c the residual bi-polaron-phonon interac- 
tion vanishes. At present we have no clear explanation 
for this free-particle like behaviour. As an example we 
plot the rescaled dispersion for different system sizes and 
diagonalization methods in the inset of Figure 3. 
(b) Luttinger liquid behaviour. The Holstein model of 
spinless fermions is defined by omitting the electron spin 
a and consequently the Coulomb interaction U in Hamil- 
tonian (1). In one dimension and at half filling, depend- 
ing on the coupling strength g, this model undergoes a 
transition from a gapless metallic phase to a Peierls dis- 
torted phase with a gap between the ground-state and 
lowest excitations. Details of this transition and the 
properties of the different phases were studied with sev- 
eral methods over the last years (see Refs. [10-14]). One 
interesting aspect is the description of the metallic phase 
in terms of an effective Luttinger model, which, accord- 
ing to the 'Luttinger liquid hypothesis' of Haldane [15], 
should be an universal picture for the low temperature 
properties of all one-dimensional metals. The two pa- 
rameters of the Luttinger model, the renormalizcd Fermi 
velocity u p and effective coupling constant K pi can be 
determined through the scaling behaviour of the ground- 
state energy E and the energy of charge excitations E±i 
with respect to the system size N: 
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In a recent work [14] we used a variational method to 
calculate eigenstates and the resulting Luttinger param- 
eters for the Holstein model at half filling. Unfortunately 
the method failed to give consistent results especially 
for K p in the anti-adiabatic regime of large frequencies 
lu 3> t where the Holstein model can be well described by 
second order perturbation theory, leading to an effective 
XXZ spin model [10] with known Luttinger parameters. 
For large frequencies the XXZ model, as well as Monte 
Carlo [12] and DMRG [13] calculations for the Holstein 
model, yield K p < 1 corresponding to a repulsive inter- 
action. If K p , starting with the value 1 for the noninter- 
acting case, reaches | with increasing coupling strength, 
the model undergoes a Kostcrlitz-Thouless transition to 
a gapped phase. In contrast, with our variational tech- 
nique we found K p > 1, corresponding to an attractive 
interaction, for all frequencies. Since it was not possible 
to calculate the required eigenstates with sufficient pre- 
cision for the various system sizes needed for finite size 
scaling, this unsatisfying situation could not be resolved 
with our 'traditional' Lanczos diagonalization procedure. 



by scaling the energies for system sizes up to N = 10. 
In the anti-adiabatic frequency range the renormalizcd 
Fermi velocity u p is drastically supressed within the 
metallic phase, while for low phonon frequencies it re- 
mains almost unchanged up to the phase transition. A 
very interesting result is the changing character of the 
interaction below u <~ t. For small frequencies the effec- 
tive fermion-fermion interaction is attractive, while it is 
repulsive for large frequencies. Possibly there is a tran- 
sition point, depending on g and u>, where the model is 
free in lowest order. Further analytical studies of this 
behaviour will be reported elsewhere. 

In conclusion, we have proposed an advanced phonon 
optimization algorithm for application in Lanczos di- 
agonalization, and demonstrated its reliability for two 
strongly interacting electron-phonon systems. 

We acknowledge valuable discussion with J. Bonca, 
H. Biittncr, E. Jeckelmann, F. Gohmann, J. Loos, 
and S.A. Trugman as well as the provision of com- 
puter resources by NIC Julich, HLR Stuttgart and LRZ 
Miinchcn. Work at Los Alamos is performed under the 
auspices of the US DOE. 
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FIG. 4.: Luttinger liquid parameters for the Holstein 
model of spinless fermions at phonon frequencies u — 
0.1, 1.0 and 10. 

We therefore reconsidered the problem with a more 
sophisticated variant of the above phonon optimization 
algorithm, using two different sets of optimized phonon 
states, one for each possible fermion occupation number 
(cf. Rcf. [5]). Together with the cut-off, this results in a 
further reduction of the Hilbert space, which is required 
for the diagonalization of larger systems. It is worth not- 
ing that this advantage is gained at the expense of a 
much more complicated fermionic Hamiltonian, since ev- 
ery hopping is connected with the projection of the actual 
phonon state onto the other basis set. Hence, for other 
models with a more difficult structure, like Jahn- Teller 
problems with two or three phonon modes per site, this 
procedure is not recommended. 

In Figure 4 we show the Luttinger parameters we found 



[1] E. K. H. Salje, A. S. Alexandrov, and W. Y. Liang, 
Polarons and Bipolarons in High Temperature Supercon- 
ductors and Related Materials, Cambridge Univ. Press, 
(Cambridge 1995). 

[2] M. Jaime, H. T. Hardner, M. R. M. B. Salamon, 
P. Dorsey, and D. Emin, Phys. Rev. Lett. 78, 951 (1997). 

[3] J. W. Bray, L. V. Interrante, I. S. Jacobs, and J. C. 
Bonner, p. 353 in Extended Linear Chain Compounds, 
ed. J. S. Miller (Plenum, New York 1985). 

[4] S.R. White, Phys. Rev. B 48, 10345 (1993). 

[5] C. Zhang, E. Jeckelmann, and S.R. White, Phys. Rev. 
Lett. 80, 2661 (1998). 

[6] B. Bauml, G. Wellein, and H. Fehske, Phys. Rev. B 58, 
3663 (1998). 

[7] A. S. Alexandrov and N. F. Mott, Rep. Prog. Phys. 57, 
1197 (1994). 

[8] J. Bonca, T. Katrasnik, S.A. Trugman, arXivxond- 
mat/0001325. 

[9] H. Fehske, H. Roder, G. Wellein, and A. Mistriotis, Phys. 
Rev. B 51, 16582 (1995). 
[10] J.E. Hirsch and E. Fradkin, Phys. Rev. B 27, 4302 
(1983). 

[11] H. Zheng, D. Feinberg, and M. Avignon, Phys. Rev. B 

39, 9405 (1989). 
[12] R.H. McKenzie, C.J. Hamer, and D.W. Murray, Phys. 

Rev. B 53, 9676 (1996). 
[13] R.J. Bursill, R.H. McKenzie, and C.J. Hamer, Phys. Rev. 

Lett. 80, 5607 (1998). 
[14] A. Weifie and H. Fehske, Phys. Rev. B 58, 13526 (1998). 
[15] F.D.M. Haldane, Phys. Rev. Lett. 45, 1358 (1980). 



4 



